Unidad Práctica 4: Visualización de Mapas¶

En este notebook trabajaremos con información geográfica. Seguiremos estudiando el data set de viajes de la encuesta origen-destino, esta vez enfocándonos en distintos patrones geográficos que puedan ayudarnos a responder preguntas específicas.

Preámbulo y Carga de Datos¶

In [1]:

import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np
import geopandas as gpd

# esto configura la calidad de la imagen. dependerá de tu resolución. el valor por omisión es 80
mpl.rcParams["figure.dpi"] = 96
# esto depende de las fuentes que tengas instaladas en el sistema.
mpl.rcParams["font.family"] = "Fira Sans Extra Condensed"

# Importar la biblioteca seaborn para la visualización de datos
import seaborn as sns

# Establecer el estilo de fondo de las gráficas como "whitegrid" en seaborn.
sns.set_style("whitegrid")

import matplotlib as mpl import matplotlib.pyplot as plt import seaborn as sns import pandas as pd import numpy as np import geopandas as gpd # esto configura la calidad de la imagen. dependerá de tu resolución. el valor por omisión es 80 mpl.rcParams["figure.dpi"] = 96 # esto depende de las fuentes que tengas instaladas en el sistema. mpl.rcParams["font.family"] = "Fira Sans Extra Condensed" # Importar la biblioteca seaborn para la visualización de datos import seaborn as sns # Establecer el estilo de fondo de las gráficas como "whitegrid" en seaborn. sns.set_style("whitegrid")

In [2]:

def decode_column(
    df,
    fname,
    col_name,
    index_col="Id",
    value_col=None,
    sep=";",
    encoding="utf-8",
    index_dtype=np.float64,
):
    """
    Agrega una columna extra al DataFrame `df` decodificando valores desde un archivo externo.

    :param df: DataFrame al que se agregará la columna extra.
    :param fname: Nombre del archivo que contiene los valores a decodificar.
    :param col_name: Nombre de la columna en el DataFrame `df` que queremos decodificar.
    :param index_col: Nombre de la columna en el archivo `fname` que contiene los índices que codifican la columna `col_name`.
    :param value_col: Nombre de la columna en el archivo `fname` que tiene los valores decodificados. Si no se proporciona, se utilizará "value" por defecto.
    :param sep: Carácter que separa los valores en el archivo `fname`.
    :param encoding: Identificación del conjunto de caracteres (character set) que utiliza el archivo. Usualmente es utf-8, si no funciona, se puede probar con iso-8859-1.
    :param index_dtype: Tipo de datos del índice en el archivo `fname`, por defecto np.float64.

    :return: Serie con los valores decodificados que se agregarán como columna extra en el DataFrame `df`.
    """
    if value_col is None:
        value_col = "value"

    # Lee el archivo `fname` que contiene los valores decodificados junto con sus índices.
    values_df = pd.read_csv(
        fname,
        sep=sep,
        index_col=index_col,
        names=[index_col, value_col],
        header=0,
        dtype={index_col: index_dtype},
        encoding=encoding,
    )

    # Extrae la columna `col_name` del DataFrame `df`.
    src_df = df.loc[:, (col_name,)]

    # Une el DataFrame original con los valores decodificados usando la columna `col_name` como índice.
    # Se obtiene una Serie con los valores decodificados que se agregarán como columna extra en el DataFrame `df`.
    return src_df.join(values_df, on=col_name)[value_col]

def decode_column( df, fname, col_name, index_col="Id", value_col=None, sep=";", encoding="utf-8", index_dtype=np.float64, ): """ Agrega una columna extra al DataFrame `df` decodificando valores desde un archivo externo. :param df: DataFrame al que se agregará la columna extra. :param fname: Nombre del archivo que contiene los valores a decodificar. :param col_name: Nombre de la columna en el DataFrame `df` que queremos decodificar. :param index_col: Nombre de la columna en el archivo `fname` que contiene los índices que codifican la columna `col_name`. :param value_col: Nombre de la columna en el archivo `fname` que tiene los valores decodificados. Si no se proporciona, se utilizará "value" por defecto. :param sep: Carácter que separa los valores en el archivo `fname`. :param encoding: Identificación del conjunto de caracteres (character set) que utiliza el archivo. Usualmente es utf-8, si no funciona, se puede probar con iso-8859-1. :param index_dtype: Tipo de datos del índice en el archivo `fname`, por defecto np.float64. :return: Serie con los valores decodificados que se agregarán como columna extra en el DataFrame `df`. """ if value_col is None: value_col = "value" # Lee el archivo `fname` que contiene los valores decodificados junto con sus índices. values_df = pd.read_csv( fname, sep=sep, index_col=index_col, names=[index_col, value_col], header=0, dtype={index_col: index_dtype}, encoding=encoding, ) # Extrae la columna `col_name` del DataFrame `df`. src_df = df.loc[:, (col_name,)] # Une el DataFrame original con los valores decodificados usando la columna `col_name` como índice. # Se obtiene una Serie con los valores decodificados que se agregarán como columna extra en el DataFrame `df`. return src_df.join(values_df, on=col_name)[value_col]

In [3]:

# leer personas
path = "https://raw.githubusercontent.com/zorzalerrante/aves/master/data/external/EOD_STGO/"
personas = pd.read_csv(
    path + "personas.csv",sep=";", decimal=",", encoding="utf-8"
    ).rename(columns={"Factor": "FactorPersona"})

# leer personas path = "https://raw.githubusercontent.com/zorzalerrante/aves/master/data/external/EOD_STGO/" personas = pd.read_csv( path + "personas.csv",sep=";", decimal=",", encoding="utf-8" ).rename(columns={"Factor": "FactorPersona"})

In [4]:

# agregar atributos 
atributos = ["Sexo","TramoIngreso","Relacion","Ocupacion"]

for col_name in atributos:
    personas[col_name] = decode_column(personas,path + f"Tablas_parametros/{col_name}.csv", col_name)

# agregar atributos atributos = ["Sexo","TramoIngreso","Relacion","Ocupacion"] for col_name in atributos: personas[col_name] = decode_column(personas,path + f"Tablas_parametros/{col_name}.csv", col_name)

In [5]:

# leer viajes
viajes  = (
        pd.read_csv(path + "viajes.csv", sep=";", decimal=",")
        .join(
            pd.read_csv(path + "ViajesDifusion.csv", sep=";", index_col="Viaje"),
            on="Viaje",
        )
        .join(
            pd.read_csv(path + "DistanciaViaje.csv", sep=";", index_col="Viaje"),
            on="Viaje",
        )
    )

# leer viajes viajes = ( pd.read_csv(path + "viajes.csv", sep=";", decimal=",") .join( pd.read_csv(path + "ViajesDifusion.csv", sep=";", index_col="Viaje"), on="Viaje", ) .join( pd.read_csv(path + "DistanciaViaje.csv", sep=";", index_col="Viaje"), on="Viaje", ) )

In [6]:

# agregar atributos 
viajes["ModoAgregado"] = decode_column(
    viajes,path + "Tablas_parametros/ModoAgregado.csv", index_col="ID",value_col="Modo",col_name = "ModoAgregado") 

viajes["ModoDifusion"] = decode_column(
    viajes,path + "Tablas_parametros/ModoDifusion.csv", encoding="latin-1",index_col="ID",col_name = "ModoDifusion") 

viajes["SectorOrigen"] = decode_column(
    viajes,path + "Tablas_parametros/Sector.csv", col_name="SectorOrigen",index_col="Sector",value_col="Nombre",sep=";")

viajes["SectorDestino"] = decode_column(
    viajes,path + "Tablas_parametros/Sector.csv", col_name="SectorDestino",index_col="Sector",value_col="Nombre",sep=";")

viajes["Proposito"] = decode_column(
    viajes,path + "Tablas_parametros/Proposito.csv",  col_name="Proposito") 

viajes["ComunaOrigen"] = decode_column(
    viajes,path + "Tablas_parametros/Comunas.csv","ComunaOrigen",value_col="Comuna",sep=",")

viajes["ComunaDestino"] = decode_column(
    viajes,path + "Tablas_parametros/Comunas.csv", "ComunaDestino",value_col="Comuna",sep=",")

viajes["Periodo"] = decode_column(
    viajes,path + "Tablas_parametros/Periodo.csv","Periodo",sep=";",value_col="Periodos")

# agregar atributos viajes["ModoAgregado"] = decode_column( viajes,path + "Tablas_parametros/ModoAgregado.csv", index_col="ID",value_col="Modo",col_name = "ModoAgregado") viajes["ModoDifusion"] = decode_column( viajes,path + "Tablas_parametros/ModoDifusion.csv", encoding="latin-1",index_col="ID",col_name = "ModoDifusion") viajes["SectorOrigen"] = decode_column( viajes,path + "Tablas_parametros/Sector.csv", col_name="SectorOrigen",index_col="Sector",value_col="Nombre",sep=";") viajes["SectorDestino"] = decode_column( viajes,path + "Tablas_parametros/Sector.csv", col_name="SectorDestino",index_col="Sector",value_col="Nombre",sep=";") viajes["Proposito"] = decode_column( viajes,path + "Tablas_parametros/Proposito.csv", col_name="Proposito") viajes["ComunaOrigen"] = decode_column( viajes,path + "Tablas_parametros/Comunas.csv","ComunaOrigen",value_col="Comuna",sep=",") viajes["ComunaDestino"] = decode_column( viajes,path + "Tablas_parametros/Comunas.csv", "ComunaDestino",value_col="Comuna",sep=",") viajes["Periodo"] = decode_column( viajes,path + "Tablas_parametros/Periodo.csv","Periodo",sep=";",value_col="Periodos")

In [7]:

# correciones
viajes = viajes[pd.notnull(viajes["HoraIni"])]
viajes = viajes[viajes["Imputada"] == 0].copy()
viajes["HoraIni"] = pd.to_timedelta(viajes["HoraIni"] + ":00")

# correciones viajes = viajes[pd.notnull(viajes["HoraIni"])] viajes = viajes[viajes["Imputada"] == 0].copy() viajes["HoraIni"] = pd.to_timedelta(viajes["HoraIni"] + ":00")

In [8]:

# Juntar Datos
viajes_persona = viajes.merge(personas)

# Juntar Datos viajes_persona = viajes.merge(personas)

In [9]:

viajes_persona["PesoLaboral"] = (
    viajes_persona["FactorLaboralNormal"] * viajes_persona["Factor_LaboralNormal"]
)

viajes_persona = viajes_persona[pd.notnull(viajes_persona["PesoLaboral"])]

print(
    "{} viajes expandidos a {}".format(
        len(viajes_persona), int(viajes_persona["PesoLaboral"].sum())
    )
)

viajes_persona["PesoLaboral"] = ( viajes_persona["FactorLaboralNormal"] * viajes_persona["Factor_LaboralNormal"] ) viajes_persona = viajes_persona[pd.notnull(viajes_persona["PesoLaboral"])] print( "{} viajes expandidos a {}".format( len(viajes_persona), int(viajes_persona["PesoLaboral"].sum()) ) )

65591 viajes expandidos a 14566372

¿Dónde se concentran las personas que utilizan cada modo de transporte en la ciudad?¶

Es de interés saber dónde viven las personas que utilizan cada modo de transporte. Eso permite informar la planificación de nuevas redes de transporte y la gestión de las redes actuales.

Esta pregunta tiene una componente geográfica: el dónde. Si tuviésemos un enfoque de tablas, podríamos utilizar el atributo Comuna o Sector y resolver tareas del tipo parte-de-un-todo para conocer la distribución del uso de modo de transporte a través de las comunas o sectores.

Sin embargo, existe variabilidad dentro de les habitantes de cada comuna. Además, también es de interés saber si dos comunas/sectores/barrios/etc. cercanos tienen comportamiento similar, y si no, entender el contexto urbano que causa la diferencia.

Ante esa necesidad, un mapa es inmejorable. Pero, ¿cómo construir y configurar el mapa? Exploremos eso ahora.

Primero, veamos cómo luces las coordenadas que vienen dentro de los datos:

In [10]:

viajes_persona[
    ["OrigenCoordX", "OrigenCoordY", "DestinoCoordX", "DestinoCoordY"]
].head()

viajes_persona[ ["OrigenCoordX", "OrigenCoordY", "DestinoCoordX", "DestinoCoordY"] ].head()

Out[10]:

	OrigenCoordX	OrigenCoordY	DestinoCoordX	DestinoCoordY
0	335208.7188	6288387.0	338812.3125	6292391.0
1	338812.2813	6292391.0	335208.7188	6288387.0
2	338536.4375	6291928.0	354267.3438	6302297.0
3	354267.3438	6302297.0	338536.4375	6291928.0
4	338536.4375	6291928.0	350841.6563	6297212.0

Si bien las coordenadas están en formato numérico, no las tenemos estructuradas en un GeoDataFrame. Para ello, usaremos la función to_point_geodataframe para darle contexto geográfico a los datos. Crearemos dos estructuras, una para los orígenes de los viajes y otra para los destinos. Del notebook anterior ya conocemos el sistema de coordenadas:

In [11]:

import shapely

def to_point_geodataframe(df, longitude, latitude, crs="epsg:4326", drop=False):
    gdf = gpd.GeoDataFrame(
        df,
        geometry=df.apply(
            lambda x: shapely.geometry.Point((x[longitude], x[latitude])), axis=1
        ),
        crs=crs,
    )

    if drop:
        gdf = gdf.drop([longitude, latitude], axis=1)

    return gdf

import shapely def to_point_geodataframe(df, longitude, latitude, crs="epsg:4326", drop=False): gdf = gpd.GeoDataFrame( df, geometry=df.apply( lambda x: shapely.geometry.Point((x[longitude], x[latitude])), axis=1 ), crs=crs, ) if drop: gdf = gdf.drop([longitude, latitude], axis=1) return gdf

In [12]:

origenes_viajes = to_point_geodataframe(
    viajes_persona, "OrigenCoordX", "OrigenCoordY", crs="epsg:32719"
)
destinos_viajes = to_point_geodataframe(
    viajes_persona, "DestinoCoordX", "DestinoCoordY", crs="epsg:32719"
)

origenes_viajes = to_point_geodataframe( viajes_persona, "OrigenCoordX", "OrigenCoordY", crs="epsg:32719" ) destinos_viajes = to_point_geodataframe( viajes_persona, "DestinoCoordX", "DestinoCoordY", crs="epsg:32719" )

Cargamos el archivo GeoJSON que creamos en el notebook 03-python-mapas-preliminario.ipynb:

In [13]:

path = "https://raw.githubusercontent.com/zorzalerrante/aves/master/data/processed/scl_zonas_urbanas.json"

zones = (
    gpd.read_file(path)
    .set_index("ID")
    .to_crs(origenes_viajes.crs)
)
zones.head()

path = "https://raw.githubusercontent.com/zorzalerrante/aves/master/data/processed/scl_zonas_urbanas.json" zones = ( gpd.read_file(path) .set_index("ID") .to_crs(origenes_viajes.crs) ) zones.head()

Out[13]:

	AREA	Zona	Com	Comuna	REGION	NOM_REGION	PROVINCIA	NOM_PROVIN	NOM_COMUNA	URBANO	TIPO	NOM_CATEG	SHAPE_Leng	SHAPE_Area	area_m2	geometry
ID
103	837.7500	103.0	13105	El Bosque	13	REGIÓN METROPOLITANA DE SANTIAGO	131	SANTIAGO	EL BOSQUE	EL BOSQUE	CAPITAL COMUNAL	CIUDAD	0.152123	0.001391	4.575649e+05	POLYGON ((346371.218 6285956.953, 346349.534 6...
104	998.8125	104.0	13105	El Bosque	13	REGIÓN METROPOLITANA DE SANTIAGO	131	SANTIAGO	EL BOSQUE	EL BOSQUE	CAPITAL COMUNAL	CIUDAD	0.152123	0.001391	7.725462e+05	POLYGON ((344415.434 6285878.048, 344548.462 6...
106	286.2500	106.0	13105	El Bosque	13	REGIÓN METROPOLITANA DE SANTIAGO	131	SANTIAGO	EL BOSQUE	EL BOSQUE	CAPITAL COMUNAL	CIUDAD	0.152123	0.001391	2.693838e+06	POLYGON ((342466.924 6284938.448, 342467.002 6...
115	857.4375	115.0	13105	El Bosque	13	REGIÓN METROPOLITANA DE SANTIAGO	131	SANTIAGO	EL BOSQUE	EL BOSQUE	CAPITAL COMUNAL	CIUDAD	0.152123	0.001391	7.534193e+05	POLYGON ((344391.421 6285910.048, 344365.314 6...
116	853.9375	116.0	13105	El Bosque	13	REGIÓN METROPOLITANA DE SANTIAGO	131	SANTIAGO	EL BOSQUE	EL BOSQUE	CAPITAL COMUNAL	CIUDAD	0.152123	0.001391	7.184305e+05	POLYGON ((345371.399 6286207.024, 345608.811 6...

Ahora nos aseguramos de tener orígenes y destinos que solo se dan dentro del contexto urbano determinado por el GeoDataFrame de zonas:

In [14]:

def clip_point_geodataframe(geodf, bounding_box, buffer=0):
    bounds = shapely.geometry.box(*bounding_box).buffer(buffer)
    return geodf[geodf.within(bounds)]

def clip_point_geodataframe(geodf, bounding_box, buffer=0): bounds = shapely.geometry.box(*bounding_box).buffer(buffer) return geodf[geodf.within(bounds)]

In [15]:

origenes_viajes = origenes_viajes[
    origenes_viajes["Viaje"].isin(destinos_viajes["Viaje"])
]
origenes_viajes = clip_point_geodataframe(origenes_viajes, zones.total_bounds)
destinos_viajes = destinos_viajes[
    destinos_viajes["Viaje"].isin(origenes_viajes["Viaje"])
]

origenes_viajes = origenes_viajes[ origenes_viajes["Viaje"].isin(destinos_viajes["Viaje"]) ] origenes_viajes = clip_point_geodataframe(origenes_viajes, zones.total_bounds) destinos_viajes = destinos_viajes[ destinos_viajes["Viaje"].isin(origenes_viajes["Viaje"]) ]

Grafiquemos los datos. Utilizaremos la función figure_from_geodataframe de aves que configura los ejes de manera automática, y el método plot de geopandas puede graficar directamente los puntos de origen/destino:

In [16]:

def figure_from_geodataframe(
    geodf,
    height=5,
    bbox=None,
    remove_axes=True,
    set_limits=True,
    basemap=None,
    basemap_interpolation="hanning",
):
    if bbox is None:
        bbox = geodf.total_bounds

    aspect = (bbox[2] - bbox[0]) / (bbox[3] - bbox[1])
    fig, ax = plt.subplots(figsize=(height * aspect, height))

    if set_limits:
        ax.set_xlim([bbox[0], bbox[2]])
        ax.set_ylim([bbox[1], bbox[3]])

        # code from geopandas
        if geodf.crs and geodf.crs.is_geographic:
            bounds = geodf.total_bounds
            y_coord = np.mean([bounds[1], bounds[3]])
            ax.set_aspect(1 / np.cos(y_coord * np.pi / 180))
            # formula ported from R package sp
            # https://github.com/edzer/sp/blob/master/R/mapasp.R
        else:
            ax.set_aspect("equal")

    if remove_axes:
        ax.set_axis_off()

    if basemap is not None:
        cx.add_basemap(
            ax,
            crs=geodf.crs.to_string(),
            source=basemap,
            interpolation=basemap_interpolation,
            zorder=0,
        )

    return fig, ax

def figure_from_geodataframe( geodf, height=5, bbox=None, remove_axes=True, set_limits=True, basemap=None, basemap_interpolation="hanning", ): if bbox is None: bbox = geodf.total_bounds aspect = (bbox[2] - bbox[0]) / (bbox[3] - bbox[1]) fig, ax = plt.subplots(figsize=(height * aspect, height)) if set_limits: ax.set_xlim([bbox[0], bbox[2]]) ax.set_ylim([bbox[1], bbox[3]]) # code from geopandas if geodf.crs and geodf.crs.is_geographic: bounds = geodf.total_bounds y_coord = np.mean([bounds[1], bounds[3]]) ax.set_aspect(1 / np.cos(y_coord * np.pi / 180)) # formula ported from R package sp # https://github.com/edzer/sp/blob/master/R/mapasp.R else: ax.set_aspect("equal") if remove_axes: ax.set_axis_off() if basemap is not None: cx.add_basemap( ax, crs=geodf.crs.to_string(), source=basemap, interpolation=basemap_interpolation, zorder=0, ) return fig, ax

In [17]:

fig, ax = figure_from_geodataframe(zones, height=9, remove_axes=True)

zones.plot(ax=ax, color="#efefef", edgecolor="#abacab", linewidth=1)
origenes_viajes.plot(ax=ax, markersize=1, marker=".", color="blue", alpha=0.05)
destinos_viajes.plot(ax=ax, markersize=1, marker=".", color="red", alpha=0.05)

plt.show()

fig, ax = figure_from_geodataframe(zones, height=9, remove_axes=True) zones.plot(ax=ax, color="#efefef", edgecolor="#abacab", linewidth=1) origenes_viajes.plot(ax=ax, markersize=1, marker=".", color="blue", alpha=0.05) destinos_viajes.plot(ax=ax, markersize=1, marker=".", color="red", alpha=0.05) plt.show()

Nos damos cuenta que hay muchos datos. Tenemos que:

• Filtrarlos, es decir, decidir cuáles son relevantes para la tarea.
• Visualizarlos utilizando la técnica adecuada.

Para filtrarlos, fíjemonos en cuatro modos de transporte: público (Bip!), privado (Auto) y activo (Caminata y Bicicleta). Y por ahora veamos los viajes al trabajo. Para visualizarlos, en aves tenemos una función dot_map que hace un mapa de puntos y que nos permitirá configurar algunos aspectos del gráfico considerando los aspectos que vimos en la clase teórica.

El gráfico se puede hacer así:

In [18]:

def bubble_map(
    ax,
    geodf: gpd.GeoDataFrame,
    size,
    scale=1,
    palette=None,
    color=None,
    add_legend=True,
    edgecolor="white",
    alpha=1.0,
    label=None,
    **kwargs
):
    marker = "o"

    if size is not None:
        if type(size) == str:
            marker_size = geodf[size]
        else:
            marker_size = float(size)
    else:
        marker_size = 1

    return geodf.plot(
        ax=ax,
        marker=marker,
        markersize=marker_size * scale,
        edgecolor=edgecolor,
        alpha=alpha,
        facecolor=color,
        legend=add_legend,
        label=label,
        **kwargs
    )


def dot_map(
    ax,
    geodf: gpd.GeoDataFrame,
    size=10,
    palette=None,
    add_legend=True,
    label=None,
    **kwargs
):
    return bubble_map(
        ax,
        geodf,
        size=float(size),
        palette=palette,
        add_legend=add_legend,
        edgecolor="none",
        label=label,
        **kwargs
    )

def bubble_map( ax, geodf: gpd.GeoDataFrame, size, scale=1, palette=None, color=None, add_legend=True, edgecolor="white", alpha=1.0, label=None, **kwargs ): marker = "o" if size is not None: if type(size) == str: marker_size = geodf[size] else: marker_size = float(size) else: marker_size = 1 return geodf.plot( ax=ax, marker=marker, markersize=marker_size * scale, edgecolor=edgecolor, alpha=alpha, facecolor=color, legend=add_legend, label=label, **kwargs ) def dot_map( ax, geodf: gpd.GeoDataFrame, size=10, palette=None, add_legend=True, label=None, **kwargs ): return bubble_map( ax, geodf, size=float(size), palette=palette, add_legend=add_legend, edgecolor="none", label=label, **kwargs )

In [19]:

fig, ax = figure_from_geodataframe(zones, height=6, remove_axes=True)

zones.plot(ax=ax, color="#efefef", edgecolor="#abacab", linewidth=1)

for idx, group in origenes_viajes[(origenes_viajes.Proposito == "Al trabajo")].groupby(
    "ModoDifusion"
):
    if idx in ["Bip!", "Auto", "Caminata", "Bicicleta"]:
        dot_map(ax, group, size=10, label=idx)

ax.legend()
plt.show()

fig, ax = figure_from_geodataframe(zones, height=6, remove_axes=True) zones.plot(ax=ax, color="#efefef", edgecolor="#abacab", linewidth=1) for idx, group in origenes_viajes[(origenes_viajes.Proposito == "Al trabajo")].groupby( "ModoDifusion" ): if idx in ["Bip!", "Auto", "Caminata", "Bicicleta"]: dot_map(ax, group, size=10, label=idx) ax.legend() plt.show()

Notamos que cada categoría utiliza un tono de color distinto. Sin embargo, cuesta ver una posible distribución geográfica debido a la oclusión y aglomeramiento de los viajes.

Una posible solución es darle un tamaño a cada punto que sea proporcional al peso de los viajes. En tal caso, en vez de dot_map utilizaremos bubble_map:

In [20]:

fig, ax = figure_from_geodataframe(zones, height=9, remove_axes=True)

zones.plot(ax=ax, color="#efefef", edgecolor="white", linewidth=1)

for idx, group in origenes_viajes[(origenes_viajes.Proposito == "Al trabajo")].groupby(
    "ModoDifusion"
):
    if idx in ["Bip!", "Auto", "Caminata", "Bicicleta"]:
        bubble_map(ax, group, size="PesoLaboral", scale=0.15, label=idx)

ax.legend()

fig.tight_layout()
plt.show()

fig, ax = figure_from_geodataframe(zones, height=9, remove_axes=True) zones.plot(ax=ax, color="#efefef", edgecolor="white", linewidth=1) for idx, group in origenes_viajes[(origenes_viajes.Proposito == "Al trabajo")].groupby( "ModoDifusion" ): if idx in ["Bip!", "Auto", "Caminata", "Bicicleta"]: bubble_map(ax, group, size="PesoLaboral", scale=0.15, label=idx) ax.legend() fig.tight_layout() plt.show()

La idea sonaba bien, sin embargo, el gráfico no nos permite responder la pregunta (aunque es un lindo gráfico para compartir y contemplar). Esto se debe a varios factores, incluyendo la interferencia entre los canales utilizados en la codificación.

También la implementación del gráfico es compleja. Requiere ciclos a través de los resultados de una operación groupby.

Utilizaremos la grilla GeoFacetGrid de aves (inspirada en FacetGrid de seaborn) para hacer el mismo gráfico pero utilizando conceptos vistos en clase.

In [21]:

from matplotlib_scalebar.scalebar import (  
    ScaleBar,
)

def north_arrow(
    ax,
    x=0.98,
    y=0.06,
    arrow_length=0.04,
    text="N",
    font_name=None,
    font_size=None,
    color="#000000",
    arrow_color="#000000",
    arrow_width=3,
    arrow_headwidth=7,
):
    ax.annotate(
        text,
        xy=(x, y),
        xytext=(x, y - arrow_length),
        arrowprops=dict(
            facecolor=arrow_color, width=arrow_width, headwidth=arrow_headwidth
        ),
        ha="center",
        va="center",
        fontsize=font_size,
        fontname=font_name,
        color=color,
        xycoords=ax.transAxes,
    )

def geographical_scale(ax, location="lower left"):
    ax.add_artist(ScaleBar(1, location="lower left"))

from matplotlib_scalebar.scalebar import ( ScaleBar, ) def north_arrow( ax, x=0.98, y=0.06, arrow_length=0.04, text="N", font_name=None, font_size=None, color="#000000", arrow_color="#000000", arrow_width=3, arrow_headwidth=7, ): ax.annotate( text, xy=(x, y), xytext=(x, y - arrow_length), arrowprops=dict( facecolor=arrow_color, width=arrow_width, headwidth=arrow_headwidth ), ha="center", va="center", fontsize=font_size, fontname=font_name, color=color, xycoords=ax.transAxes, ) def geographical_scale(ax, location="lower left"): ax.add_artist(ScaleBar(1, location="lower left"))

In [22]:

from seaborn.axisgrid import FacetGrid, Grid
from matplotlib import colorbar


class GeoFacetGrid(FacetGrid):
    def __init__(self, geodataframe: gpd.GeoDataFrame, *args, **kwargs):
        geocontext = kwargs.pop("context", None)

        if geocontext is None:
            geocontext = geodataframe

        self.geocontext = geocontext

        self.bounds = geocontext.total_bounds
        self.aspect = (self.bounds[2] - self.bounds[0]) / (
            self.bounds[3] - self.bounds[1]
        )

        kwargs["aspect"] = self.aspect
        kwargs["xlim"] = (self.bounds[0], self.bounds[2])
        kwargs["ylim"] = (self.bounds[1], self.bounds[3])

        super().__init__(geodataframe, *args, **kwargs)

        for ax in self.axes.flatten():
            if kwargs.get("remove_axes", True):
                ax.set_axis_off()

            aspect = kwargs.get("aspect", "auto")
            # code from geopandas
            if aspect == "auto":
                if geodataframe.crs and geodataframe.crs.is_geographic:
                    bounds = geodataframe.total_bounds
                    y_coord = np.mean([bounds[1], bounds[3]])
                    ax.set_aspect(1 / np.cos(y_coord * np.pi / 180))
                    # formula ported from R package sp
                    # https://github.com/edzer/sp/blob/master/R/mapasp.R
                else:
                    ax.set_aspect("equal")
            elif aspect is not None:
                ax.set_aspect(aspect)

        self.zorder = 0

    def add_layer(self, func_or_data, *args, **kwargs):
        if isinstance(func_or_data, gpd.GeoDataFrame):
            # a direct geography
            for ax in self.axes.flatten():
                func_or_data.plot(*args, ax=ax, zorder=self.zorder, **kwargs)
        else:
            plot = lambda *a, **kw: func_or_data(
                plt.gca(), kw.pop("data"), *a, zorder=self.zorder, **kw
            )
            self.map_dataframe(plot, *args, **kwargs)

        self.zorder += 1

    def add_basemap(
        self, file_path, interpolation="hanning", reset_extent=False, **kwargs
    ):
        for ax in self.axes.flatten():
            cx.add_basemap(
                ax,
                crs=self.geocontext.crs.to_string(),
                source=file_path,
                interpolation=interpolation,
                zorder=self.zorder,
                reset_extent=reset_extent,
                **kwargs
            )

            if not reset_extent:
                ax.set_xlim((self.bounds[0], self.bounds[2]))
                ax.set_ylim((self.bounds[1], self.bounds[3]))

        self.zorder += 1

    def add_map_elements(
        self,
        all_axes=False,
        scale=True,
        arrow=True,
        scale_args={},
        arrow_args={},
    ):
        for ax in self.axes.flatten():
            if arrow:
                north_arrow(ax, **arrow_args)

            if scale:
                geographical_scale(ax, **scale_args)
            if not all_axes:
                break

    def add_global_colorbar(self, palette, k, title=None, title_args={}, **kwargs):

        orientation = kwargs.get("orientation", "horizontal")
        if orientation == "horizontal":
            cax = self.fig.add_axes([0.25, -0.012, 0.5, 0.01])
        elif orientation == "vertical":
            cax = self.fig.add_axes([1.01, 0.25, 0.01, 0.5])
        else:
            raise ValueError("unsupported colorbar orientation")

        if title:
            cax.set_title(title, **title_args)

        cb = colorbar.ColorbarBase(
            cax, cmap=colormap_from_palette(palette, n_colors=k), **kwargs
        )

        cax.set_axis_off()

        return cax, cb

    def set_title(self, title, **kwargs):
        self.fig.suptitle(title, **kwargs)

from seaborn.axisgrid import FacetGrid, Grid from matplotlib import colorbar class GeoFacetGrid(FacetGrid): def __init__(self, geodataframe: gpd.GeoDataFrame, *args, **kwargs): geocontext = kwargs.pop("context", None) if geocontext is None: geocontext = geodataframe self.geocontext = geocontext self.bounds = geocontext.total_bounds self.aspect = (self.bounds[2] - self.bounds[0]) / ( self.bounds[3] - self.bounds[1] ) kwargs["aspect"] = self.aspect kwargs["xlim"] = (self.bounds[0], self.bounds[2]) kwargs["ylim"] = (self.bounds[1], self.bounds[3]) super().__init__(geodataframe, *args, **kwargs) for ax in self.axes.flatten(): if kwargs.get("remove_axes", True): ax.set_axis_off() aspect = kwargs.get("aspect", "auto") # code from geopandas if aspect == "auto": if geodataframe.crs and geodataframe.crs.is_geographic: bounds = geodataframe.total_bounds y_coord = np.mean([bounds[1], bounds[3]]) ax.set_aspect(1 / np.cos(y_coord * np.pi / 180)) # formula ported from R package sp # https://github.com/edzer/sp/blob/master/R/mapasp.R else: ax.set_aspect("equal") elif aspect is not None: ax.set_aspect(aspect) self.zorder = 0 def add_layer(self, func_or_data, *args, **kwargs): if isinstance(func_or_data, gpd.GeoDataFrame): # a direct geography for ax in self.axes.flatten(): func_or_data.plot(*args, ax=ax, zorder=self.zorder, **kwargs) else: plot = lambda *a, **kw: func_or_data( plt.gca(), kw.pop("data"), *a, zorder=self.zorder, **kw ) self.map_dataframe(plot, *args, **kwargs) self.zorder += 1 def add_basemap( self, file_path, interpolation="hanning", reset_extent=False, **kwargs ): for ax in self.axes.flatten(): cx.add_basemap( ax, crs=self.geocontext.crs.to_string(), source=file_path, interpolation=interpolation, zorder=self.zorder, reset_extent=reset_extent, **kwargs ) if not reset_extent: ax.set_xlim((self.bounds[0], self.bounds[2])) ax.set_ylim((self.bounds[1], self.bounds[3])) self.zorder += 1 def add_map_elements( self, all_axes=False, scale=True, arrow=True, scale_args={}, arrow_args={}, ): for ax in self.axes.flatten(): if arrow: north_arrow(ax, **arrow_args) if scale: geographical_scale(ax, **scale_args) if not all_axes: break def add_global_colorbar(self, palette, k, title=None, title_args={}, **kwargs): orientation = kwargs.get("orientation", "horizontal") if orientation == "horizontal": cax = self.fig.add_axes([0.25, -0.012, 0.5, 0.01]) elif orientation == "vertical": cax = self.fig.add_axes([1.01, 0.25, 0.01, 0.5]) else: raise ValueError("unsupported colorbar orientation") if title: cax.set_title(title, **title_args) cb = colorbar.ColorbarBase( cax, cmap=colormap_from_palette(palette, n_colors=k), **kwargs ) cax.set_axis_off() return cax, cb def set_title(self, title, **kwargs): self.fig.suptitle(title, **kwargs)

In [23]:

#from aves.visualization.figures import GeoFacetGrid

grid = GeoFacetGrid(
    # los datos
    origenes_viajes,
    # el contexto geográfico. es opcional, se utiliza para configurar el gráfico
    context=zones,
    # la variable que mapearemos a las filas del gráfico
    row="Proposito",
    # en este caso, solo una fila
    row_order=["Al trabajo"],
    # el canal de codificación tono (hue) expresará la columna categórica ModoDifusion
    hue="ModoDifusion",
    # solo veremos estas categorías
    hue_order=["Auto", "Bicicleta", "Bip!", "Caminata"],
    # la altura del gráfico
    height=9,
)

# visualizamos el contexto
grid.add_layer(zones, color="#efefef", edgecolor="white", linewidth=1)

# agregamos los bubble_map correspondientes
grid.add_layer(bubble_map, size="PesoLaboral", scale=0.15)

# agregamos la leyenda
grid.add_legend()

grid.fig.tight_layout()

#from aves.visualization.figures import GeoFacetGrid grid = GeoFacetGrid( # los datos origenes_viajes, # el contexto geográfico. es opcional, se utiliza para configurar el gráfico context=zones, # la variable que mapearemos a las filas del gráfico row="Proposito", # en este caso, solo una fila row_order=["Al trabajo"], # el canal de codificación tono (hue) expresará la columna categórica ModoDifusion hue="ModoDifusion", # solo veremos estas categorías hue_order=["Auto", "Bicicleta", "Bip!", "Caminata"], # la altura del gráfico height=9, ) # visualizamos el contexto grid.add_layer(zones, color="#efefef", edgecolor="white", linewidth=1) # agregamos los bubble_map correspondientes grid.add_layer(bubble_map, size="PesoLaboral", scale=0.15) # agregamos la leyenda grid.add_legend() grid.fig.tight_layout()

La implementación del gráfico es distinta, sin embargo, los problemas conceptuales se mantienen. Ahora bien, al utilizar la grilla, podemos utilizar sus opciones para reducir su complejidad a través de facetas.

Por ejemplo, podemos dividir el gráfico en columnas, donde cada columna muestra un modo de transporte diferente:

In [24]:

grid = GeoFacetGrid(
    origenes_viajes,
    context=zones,
    row="Proposito",
    col="ModoDifusion",
    row_order=["Al trabajo", "De salud"],
    col_order=["Auto", "Bip!", "Caminata", "Bicicleta"],
    height=9,
    hue="ModoDifusion",
    palette="plasma",
)

grid.add_layer(zones, color="#efefef", edgecolor="white", linewidth=1)
grid.add_layer(bubble_map, size="PesoLaboral", scale=0.15, edgecolor="black", alpha=0.5)

grid = GeoFacetGrid( origenes_viajes, context=zones, row="Proposito", col="ModoDifusion", row_order=["Al trabajo", "De salud"], col_order=["Auto", "Bip!", "Caminata", "Bicicleta"], height=9, hue="ModoDifusion", palette="plasma", ) grid.add_layer(zones, color="#efefef", edgecolor="white", linewidth=1) grid.add_layer(bubble_map, size="PesoLaboral", scale=0.15, edgecolor="black", alpha=0.5)

¡Mucho mejor! Sin embargo, la sobreposición de las burbujas sigue siendo un problema. No logramos ver la distribución debido a la oclusión de los puntos.

La solución es usar un heat_map. Así se ve un heatmap sobre todos los viajes al trabajo. Solo debemos reemplazar el uso de bubble_map por heat_map (y usar los parámetros correspondientes):

In [25]:

import KDEpy
import matplotlib.colors as colors

def kde_from_points(
    geodf,
    kernel="gaussian",
    norm=2,
    bandwidth=1e-2,
    grid_points=2 ** 9,
    weight_column=None,
):
    # La variable grid_points define la cantidad de puntos en el espacio en el que se estimará la densidad
    # hacemos una lista con las coordenadas de los viajes
    point_coords = np.vstack([geodf.geometry.x, geodf.geometry.y]).T
    # instanciamos la Fast-Fourier Transform Kernel Density Estimation
    kde = KDEpy.FFTKDE(bw=bandwidth, norm=norm, kernel=kernel)
    weights = None if weight_column is None else geodf[weight_column].values
    grid, points = kde.fit(point_coords, weights=weights).evaluate(grid_points)
    x, y = np.unique(grid[:, 0]), np.unique(grid[:, 1])
    z = points.reshape(grid_points, grid_points).T
    return x, y, z

import KDEpy import matplotlib.colors as colors def kde_from_points( geodf, kernel="gaussian", norm=2, bandwidth=1e-2, grid_points=2 ** 9, weight_column=None, ): # La variable grid_points define la cantidad de puntos en el espacio en el que se estimará la densidad # hacemos una lista con las coordenadas de los viajes point_coords = np.vstack([geodf.geometry.x, geodf.geometry.y]).T # instanciamos la Fast-Fourier Transform Kernel Density Estimation kde = KDEpy.FFTKDE(bw=bandwidth, norm=norm, kernel=kernel) weights = None if weight_column is None else geodf[weight_column].values grid, points = kde.fit(point_coords, weights=weights).evaluate(grid_points) x, y = np.unique(grid[:, 0]), np.unique(grid[:, 1]) z = points.reshape(grid_points, grid_points).T return x, y, z

In [26]:

def colormap_from_palette(palette_name, n_colors=10):
    return colors.ListedColormap(sns.color_palette(palette_name, n_colors=n_colors))

def colormap_from_palette(palette_name, n_colors=10): return colors.ListedColormap(sns.color_palette(palette_name, n_colors=n_colors))

In [27]:

def heat_map(
    ax,
    geodf,
    weight=None,
    low_threshold=0,
    max_threshold=1.0,
    n_levels=5,
    alpha=1.0,
    palette="magma",
    kernel="cosine",
    norm=2,
    bandwidth=1e-2,
    grid_points=2 ** 6,
    return_heat=False,
    cbar_label=None,
    cbar_width=2.4,
    cbar_height=0.15,
    cbar_location="upper left",
    cbar_orientation="horizontal",
    cbar_pad=0.05,
    cbar_bbox_to_anchor=(0.0, 0.0, 1.0, 1.0),
    cbar_bbox_transform=None,
    legend_type="none",
    **kwargs
):
    heat = kde_from_points(
        geodf,
        kernel=kernel,
        norm=norm,
        bandwidth=bandwidth,
        grid_points=grid_points,
        weight_column=weight,
    )

    norm_heat = heat[2] / heat[2].max()
    cmap = colormap_from_palette(palette, n_colors=n_levels)

    levels = np.linspace(low_threshold, max_threshold, n_levels)

    # TODO: this should be factorized into an utility function
    if legend_type == "colorbar":
        # add_ranged_color_legend(ax)
        if cbar_bbox_transform is None:
            cbar_bbox_transform = ax.transAxes

        if cbar_location != "out":
            cbar_ax = inset_axes(
                ax,
                width=cbar_width,
                height=cbar_height,
                loc=cbar_location,
                bbox_to_anchor=cbar_bbox_to_anchor,
                bbox_transform=cbar_bbox_transform,
                borderpad=0,
            )
        else:
            divider = make_axes_locatable(ax)
            cbar_main = divider.append_axes(
                "bottom" if cbar_orientation == "horizontal" else "right",
                size=cbar_height if cbar_orientation == "horizontal" else cbar_width,
                pad=cbar_pad,
            )
            cbar_main.set_axis_off()
            cbar_ax = inset_axes(
                cbar_main,
                width=cbar_width,
                height=cbar_height,
                loc="center",
                bbox_to_anchor=(0.0, 0.0, 1.0, 1.0),
                bbox_transform=cbar_main.transAxes,
                borderpad=0,
            )

        color_legend(cbar_ax, cmap, levels, orientation=cbar_orientation)
    else:
        cbar_ax = None

    if not return_heat:
        return (
            ax.contourf(heat[0], heat[1], norm_heat, levels, alpha=alpha, cmap=cmap),
            cbar_ax,
        )
    else:
        return (
            ax.contourf(heat[0], heat[1], norm_heat, levels, alpha=alpha, cmap=cmap),
            cbar_ax,
            heat,
        )

def heat_map( ax, geodf, weight=None, low_threshold=0, max_threshold=1.0, n_levels=5, alpha=1.0, palette="magma", kernel="cosine", norm=2, bandwidth=1e-2, grid_points=2 ** 6, return_heat=False, cbar_label=None, cbar_width=2.4, cbar_height=0.15, cbar_location="upper left", cbar_orientation="horizontal", cbar_pad=0.05, cbar_bbox_to_anchor=(0.0, 0.0, 1.0, 1.0), cbar_bbox_transform=None, legend_type="none", **kwargs ): heat = kde_from_points( geodf, kernel=kernel, norm=norm, bandwidth=bandwidth, grid_points=grid_points, weight_column=weight, ) norm_heat = heat[2] / heat[2].max() cmap = colormap_from_palette(palette, n_colors=n_levels) levels = np.linspace(low_threshold, max_threshold, n_levels) # TODO: this should be factorized into an utility function if legend_type == "colorbar": # add_ranged_color_legend(ax) if cbar_bbox_transform is None: cbar_bbox_transform = ax.transAxes if cbar_location != "out": cbar_ax = inset_axes( ax, width=cbar_width, height=cbar_height, loc=cbar_location, bbox_to_anchor=cbar_bbox_to_anchor, bbox_transform=cbar_bbox_transform, borderpad=0, ) else: divider = make_axes_locatable(ax) cbar_main = divider.append_axes( "bottom" if cbar_orientation == "horizontal" else "right", size=cbar_height if cbar_orientation == "horizontal" else cbar_width, pad=cbar_pad, ) cbar_main.set_axis_off() cbar_ax = inset_axes( cbar_main, width=cbar_width, height=cbar_height, loc="center", bbox_to_anchor=(0.0, 0.0, 1.0, 1.0), bbox_transform=cbar_main.transAxes, borderpad=0, ) color_legend(cbar_ax, cmap, levels, orientation=cbar_orientation) else: cbar_ax = None if not return_heat: return ( ax.contourf(heat[0], heat[1], norm_heat, levels, alpha=alpha, cmap=cmap), cbar_ax, ) else: return ( ax.contourf(heat[0], heat[1], norm_heat, levels, alpha=alpha, cmap=cmap), cbar_ax, heat, )

In [28]:

grid = GeoFacetGrid(
    origenes_viajes,
    context=zones,
    row="Proposito",
    col="ModoDifusion",
    row_order=["Al trabajo"],
    col_order=["Auto", "Bip!", "Caminata", "Bicicleta"],
    height=9,
    hue="ModoDifusion"
)

grid.add_layer(zones, color="#efefef", edgecolor="white", linewidth=1)

grid.add_layer(
    heat_map,
    # atributo de los datos con la importancia o peso de cada viaje
    weight="PesoLaboral",
    # cantidad de niveles/colores del mapa de calor
    n_levels=10,
    # radio de influencia de cada viaje
    bandwidth=1000,
    # valor de corte para los valores bajos del heatmap
    low_threshold=0.025,
    # transparencia
    alpha=0.75,
    # paleta de colores
    palette="inferno"
)

grid = GeoFacetGrid( origenes_viajes, context=zones, row="Proposito", col="ModoDifusion", row_order=["Al trabajo"], col_order=["Auto", "Bip!", "Caminata", "Bicicleta"], height=9, hue="ModoDifusion" ) grid.add_layer(zones, color="#efefef", edgecolor="white", linewidth=1) grid.add_layer( heat_map, # atributo de los datos con la importancia o peso de cada viaje weight="PesoLaboral", # cantidad de niveles/colores del mapa de calor n_levels=10, # radio de influencia de cada viaje bandwidth=1000, # valor de corte para los valores bajos del heatmap low_threshold=0.025, # transparencia alpha=0.75, # paleta de colores palette="inferno" )

Esta técnica de visualización responde nuestra pregunta: nos dice dónde hay concentración, y, al mismo tiempo, deduce una forma, unos límites, para esos lugares que concentran actividad.

Pero hay algo extraño en el gráfico. Solo vemos concentración en uno o dos lugares de la ciudad. Queríamos analizar de acuerdo al lugar de residencia de las personas, en un tipo de viaje importante, pero no estamos viendo mapas que nos hagan pensar en residencia. ¿Qué sucede?

La respuesta es que nos fijamos en los orígenes de los viajes al trabajo, pero sin discriminar que el viaje al trabajo fuese desde el hogar de las personas. Y, en un día laboral normal, usualmente se hacen dos viajes al trabajo. Uno desde la casa, y otro al regresar al trabajo luego de comer.

Hay soluciones para esto. La primera sería corregir ese análisis. No lo haremos puesto que no es una clase de análisis de datos. Queda como problema propuesto.

La segunda es utilizar los destinos de los viajes. En un día laboral normal podemos asumir que es común que solo se regrese a casa una vez, y que ese regreso se realice en el modo de transporte principal.

In [29]:

grid = GeoFacetGrid(
    destinos_viajes,
    context=zones,
    col="ModoDifusion",
    row="Proposito",
    row_order=["volver a casa"],
    col_order=["Auto", "Bip!", "Caminata", "Bicicleta"],
    height=9,
)

grid.add_layer(zones, color="#efefef", edgecolor="white", linewidth=1)
grid.add_layer(
    heat_map,
    weight="PesoLaboral",
    n_levels=10,
    bandwidth=1000,
    low_threshold=0.025,
    alpha=0.75,
    palette="inferno",
)

grid = GeoFacetGrid( destinos_viajes, context=zones, col="ModoDifusion", row="Proposito", row_order=["volver a casa"], col_order=["Auto", "Bip!", "Caminata", "Bicicleta"], height=9, ) grid.add_layer(zones, color="#efefef", edgecolor="white", linewidth=1) grid.add_layer( heat_map, weight="PesoLaboral", n_levels=10, bandwidth=1000, low_threshold=0.025, alpha=0.75, palette="inferno", )

Configuremos el gráfico para que sea más amigable. Haremos lo siguiente:

• Haremos un gráfico de 2x2. Para ello necesitamos filtrar por nuestra cuenta los datos, de modo que el objeto grid pueda disponer de las filas (o de las columnas) repitiendo valores.
• Agregaremos como base la imagen del territorio que descargamos en el notebook anterior.
• Utilizaremos una leyenda de colores para explicar cuál es el color asociado a la concentración.

In [30]:

import contextily as cx

grid = GeoFacetGrid(
    destinos_viajes[destinos_viajes["Proposito"] == "volver a casa"],
    context=zones,
    col="ModoDifusion",
    col_wrap=2,
    col_order=["Auto", "Bip!", "Caminata", "Bicicleta"],
    height=9,
)

grid.add_basemap("data/scl_toner_12.tif")
grid.add_layer(
    heat_map,
    weight="PesoLaboral",
    n_levels=10,
    bandwidth=1000,
    low_threshold=0.025,
    alpha=0.75,
    palette="inferno",
)
grid.add_global_colorbar(
    "inferno",
    10,
    title="Intensidad de Viajes (de menos a más)",
    orientation="horizontal",
)
grid.tight_layout()
# grid.savefig('../reports/figures/example_geofacetgrid.png', dpi=150)

import contextily as cx grid = GeoFacetGrid( destinos_viajes[destinos_viajes["Proposito"] == "volver a casa"], context=zones, col="ModoDifusion", col_wrap=2, col_order=["Auto", "Bip!", "Caminata", "Bicicleta"], height=9, ) grid.add_basemap("data/scl_toner_12.tif") grid.add_layer( heat_map, weight="PesoLaboral", n_levels=10, bandwidth=1000, low_threshold=0.025, alpha=0.75, palette="inferno", ) grid.add_global_colorbar( "inferno", 10, title="Intensidad de Viajes (de menos a más)", orientation="horizontal", ) grid.tight_layout() # grid.savefig('../reports/figures/example_geofacetgrid.png', dpi=150)

C:\Users\franc\AppData\Local\pypoetry\Cache\virtualenvs\db-connectors-RsEieBu8-py3.8\lib\site-packages\seaborn\axisgrid.py:118: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  self._figure.tight_layout(*args, **kwargs)

Out[30]:

<__main__.GeoFacetGrid at 0x18c0e988550>

Agregar el contexto territorial nos permite ver cosas como el rol segregador que juega la autopista Américo Vespucio. ¿Qué otras cosas ven ustedes?

Ejercicio propuesto:

• Consideren el análisis por grupos etáreos que hicimos la clase pasada. Definan categorías de grupos y determinen cuáles lugares de la ciudad visitan, tanto en general como por propósito de viaje. Rangos de ejemplo incluyen: menores de edad (< 18 años), adultos mayores (> 65 años), etc.

Nota: Debes tener descargado el archivo scl_toner_12.tif en tu computador local o google colab.

¿Cuán lejos queda el trabajo de acuerdo al lugar de residencia?¶

Con esta pregunta queremos entender si existe un patrón geográfico en las elecciones de residencia y trabajo de las personas.

Para responder la pregunta, primero filtramos los viajes que nos interesan:

In [31]:

viajes_trabajo = origenes_viajes[(origenes_viajes.Proposito == 'Al trabajo') &
                                (pd.notnull(origenes_viajes.PesoLaboral)) &
                                (origenes_viajes.DistEuclidiana > 0)].drop_duplicates(subset='Persona', keep='first')
                                
print(len(viajes_trabajo), viajes_trabajo.PesoLaboral.sum())

viajes_trabajo = origenes_viajes[(origenes_viajes.Proposito == 'Al trabajo') & (pd.notnull(origenes_viajes.PesoLaboral)) & (origenes_viajes.DistEuclidiana > 0)].drop_duplicates(subset='Persona', keep='first') print(len(viajes_trabajo), viajes_trabajo.PesoLaboral.sum())

11238 2380126.1986132185

Observamos que, si bien son 11228 filas en la tabla.

La columna DistEuclidiana contiene la distancia entre los puntos de origen y destino de los viajes.

No todos los viajes al trabajo se efectúan desde la casa, ni una persona hace un único viaje al trabajo durante el día. Como una manera de alivianar el problema nos quedamos con un único viaje al trabajo por persona.

Veamos la distribución de la distancia:

In [32]:

viajes_trabajo['DistEuclidiana'].describe()

viajes_trabajo['DistEuclidiana'].describe()

Out[32]:

count    1.123800e+04
mean     1.008620e+04
std      1.900365e+04
min      1.000000e+00
25%      3.628250e+03
50%      8.025000e+03
75%      1.414975e+04
max      1.258473e+06
Name: DistEuclidiana, dtype: float64

Ahora bien, el promedio que nos entregó el método pd.describe no considera la representatividad de cada viaje. Podemos utilizar la función weighted_mean para calcular el promedio ponderado:

In [33]:

def weighted_mean(df, value_column, weighs_column):
    weighted_sum = (df[value_column] * df[weighs_column]).sum()
    return weighted_sum / df[weighs_column].sum()

def weighted_mean(df, value_column, weighs_column): weighted_sum = (df[value_column] * df[weighs_column]).sum() return weighted_sum / df[weighs_column].sum()

In [34]:

viajes_trabajo['DistEuclidiana'].mean(), weighted_mean(viajes_trabajo, 'DistEuclidiana', 'PesoLaboral')

viajes_trabajo['DistEuclidiana'].mean(), weighted_mean(viajes_trabajo, 'DistEuclidiana', 'PesoLaboral')

Out[34]:

(10086.19514148425, 9549.626344374019)

En este caso, el promedio ponderado no está tan lejos del promedio sin ponderar, pero nada asegura que esa diferencia sea pequeña en toda la ciudad.

Podemos calcular la distancia promedio al trabajo por zona haciendo una operación pd.groupby:

In [35]:

distancia_zonas = (viajes_trabajo
                   .groupby(['ZonaOrigen'])
                   .apply(lambda x: weighted_mean(x, 'DistEuclidiana', 'PesoLaboral'))
                   .rename('distancia_al_trabajo')
)

distancia_zonas

distancia_zonas = (viajes_trabajo .groupby(['ZonaOrigen']) .apply(lambda x: weighted_mean(x, 'DistEuclidiana', 'PesoLaboral')) .rename('distancia_al_trabajo') ) distancia_zonas

Out[35]:

ZonaOrigen
1       6452.138966
2       5952.584070
3       4091.135703
4       4731.477408
5      17567.000000
           ...     
857    22153.295597
858    18184.397206
859    17326.810773
860    16946.779972
861    12782.225811
Name: distancia_al_trabajo, Length: 738, dtype: float64

La serie contiene, para cada zona, la distancia promedio al trabajo de la gente que vive en ella (y posiblemente otras más). Si queremos ver esta distribución podemos utilizar el método plot(kind='kde'), como hicimos en notebooks anteriores:

In [36]:

distancia_zonas.plot(kind='kde')
plt.xlim([0, distancia_zonas.max()])
plt.title('Distancia al Trabajo por Zonas')
plt.xlabel('Distancia')
plt.ylabel('Densidad (KDE)')
sns.despine()

distancia_zonas.plot(kind='kde') plt.xlim([0, distancia_zonas.max()]) plt.title('Distancia al Trabajo por Zonas') plt.xlabel('Distancia') plt.ylabel('Densidad (KDE)') sns.despine()

Ahora bien, esta distribución solamente nos permite saber propiedades estádisticas. ¡No tenemos un contexto geográfico! Por eso necesitamos el mapa, para saber si la distancia tiene relación con la ubicación de cada zona en la ciudad, con saber si zonas que están cerca entre sí tienen distancias al trabajo similares.

Como tenemos los valores de distancia por zona, un choropleth map nos permitiría mostrar la distancia asociada a cada zona al mismo tiempo que su posición geográfica.

In [37]:

zones.crs

zones.crs

Out[37]:

<Projected CRS: EPSG:32719>
Name: WGS 84 / UTM zone 19S
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Between 72°W and 66°W, southern hemisphere between 80°S and equator, onshore and offshore. Argentina. Bolivia. Brazil. Chile. Colombia. Peru.
- bounds: (-72.0, -80.0, -66.0, 0.0)
Coordinate Operation:
- name: UTM zone 19S
- method: Transverse Mercator
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

In [38]:

def color_legend(
    ax,
    color_list,
    bins=None,
    sizes=None,
    orientation="horizontal",
    remove_axes=False,
    bin_spacing="proportional",
    tick_labels=None,
):
    if bins is None:
        if type(color_list) == colors.ListedColormap:
            N = color_list.N
        else:
            N = len(color_list)
        bins = np.array(range(N))

    if sizes is not None:
        bar_width = bins[1:] - bins[0:-1]
        if orientation == "horizontal":
            ax.bar(
                bins[:-1],
                sizes,
                width=bar_width,
                align="edge",
                color=color_list,
                edgecolor=color_list,
            )
            ax.set_xticks(bins, labels=tick_labels)
        else:
            ax.barh(
                bins[:-1],
                sizes,
                height=bar_width,
                align="edge",
                color=color_list,
                edgecolor=color_list,
            )
            ax.set_yticks(bins, labels=tick_labels)
    else:
        cbar_norm = colors.BoundaryNorm(bins, len(bins) - 1)
        if type(color_list) == colors.ListedColormap:
            cmap = color_list
        else:
            cmap = colors.ListedColormap(color_list)
        cb = colorbar.ColorbarBase(
            ax,
            cmap=cmap,
            norm=cbar_norm,
            ticks=bins,
            spacing=bin_spacing,
            orientation=orientation,
        )
        if tick_labels:
            cb.set_ticklabels(tick_labels)

    sns.despine(ax=ax, top=True, bottom=True, left=True, right=True)

    if remove_axes:
        ax.set_axis_off()

    return ax

def color_legend( ax, color_list, bins=None, sizes=None, orientation="horizontal", remove_axes=False, bin_spacing="proportional", tick_labels=None, ): if bins is None: if type(color_list) == colors.ListedColormap: N = color_list.N else: N = len(color_list) bins = np.array(range(N)) if sizes is not None: bar_width = bins[1:] - bins[0:-1] if orientation == "horizontal": ax.bar( bins[:-1], sizes, width=bar_width, align="edge", color=color_list, edgecolor=color_list, ) ax.set_xticks(bins, labels=tick_labels) else: ax.barh( bins[:-1], sizes, height=bar_width, align="edge", color=color_list, edgecolor=color_list, ) ax.set_yticks(bins, labels=tick_labels) else: cbar_norm = colors.BoundaryNorm(bins, len(bins) - 1) if type(color_list) == colors.ListedColormap: cmap = color_list else: cmap = colors.ListedColormap(color_list) cb = colorbar.ColorbarBase( ax, cmap=cmap, norm=cbar_norm, ticks=bins, spacing=bin_spacing, orientation=orientation, ) if tick_labels: cb.set_ticklabels(tick_labels) sns.despine(ax=ax, top=True, bottom=True, left=True, right=True) if remove_axes: ax.set_axis_off() return ax

In [39]:

from mpl_toolkits.axes_grid1.inset_locator import inset_axes

def add_ranged_color_legend(
    ax,
    bins,
    built_palette,
    location="lower center",
    orientation="horizontal",
    label=None,
    label_size="x-small",
    title_size="small",
    title_align="left",
    width="50%",
    height="5%",
    bbox_to_anchor=(0.0, 0.0, 0.95, 0.95),
    bbox_transform=None,
    tick_labels=None,
    **kwargs,
):
    if bbox_transform is None:
        bbox_transform = ax.transAxes

    if location != "out":
        cbar_ax = inset_axes(
            ax,
            width=width,
            height=height,
            loc=location,
            bbox_to_anchor=bbox_to_anchor,
            bbox_transform=bbox_transform,
        )
    else:
        divider = make_axes_locatable(ax)
        cbar_main = divider.append_axes(
            "bottom" if orientation == "horizontal" else "right",
            size=height if orientation == "horizontal" else width,
        )
        cbar_main.set_axis_off()
        cbar_ax = inset_axes(
            cbar_main,
            width=width,
            height=height,
            loc="center",
            bbox_to_anchor=(0.0, 0.0, 1.0, 1.0),
            bbox_transform=cbar_main.transAxes,
            borderpad=0,
        )

    color_legend(
        cbar_ax,
        built_palette,
        bins,
        orientation=orientation,
        remove_axes=False,
        tick_labels=tick_labels if tick_labels is not None else [],
        **kwargs,
    )

    sns.despine(ax=cbar_ax, left=True, top=True, bottom=True, right=True)
    if label:
        cbar_ax.set_title(label, loc="left", fontsize=title_size)

    cbar_ax.tick_params(labelsize=label_size)

    return cbar_ax

from mpl_toolkits.axes_grid1.inset_locator import inset_axes def add_ranged_color_legend( ax, bins, built_palette, location="lower center", orientation="horizontal", label=None, label_size="x-small", title_size="small", title_align="left", width="50%", height="5%", bbox_to_anchor=(0.0, 0.0, 0.95, 0.95), bbox_transform=None, tick_labels=None, **kwargs, ): if bbox_transform is None: bbox_transform = ax.transAxes if location != "out": cbar_ax = inset_axes( ax, width=width, height=height, loc=location, bbox_to_anchor=bbox_to_anchor, bbox_transform=bbox_transform, ) else: divider = make_axes_locatable(ax) cbar_main = divider.append_axes( "bottom" if orientation == "horizontal" else "right", size=height if orientation == "horizontal" else width, ) cbar_main.set_axis_off() cbar_ax = inset_axes( cbar_main, width=width, height=height, loc="center", bbox_to_anchor=(0.0, 0.0, 1.0, 1.0), bbox_transform=cbar_main.transAxes, borderpad=0, ) color_legend( cbar_ax, built_palette, bins, orientation=orientation, remove_axes=False, tick_labels=tick_labels if tick_labels is not None else [], **kwargs, ) sns.despine(ax=cbar_ax, left=True, top=True, bottom=True, right=True) if label: cbar_ax.set_title(label, loc="left", fontsize=title_size) cbar_ax.tick_params(labelsize=label_size) return cbar_ax

In [40]:

from mapclassify import FisherJenks, Quantiles

def choropleth_map(
    ax,
    geodf,
    column,
    k=10,
    palette=None,
    default_divergent="RdBu_r",
    default_negative="Blues_r",
    default_positive="Reds",
    palette_type="light",
    legend="colorbar",
    edgecolor="white",
    palette_center=None,
    binning="uniform",
    bins=None,
    alpha=1.0,
    linewidth=1,
    zorder=1,
    cbar_args={},
    **kwargs,
):
    geodf = geodf[pd.notnull(geodf[column])].copy()
    min_value, max_value = geodf[column].min(), geodf[column].max()

    if binning in ("fisher_jenks", "quantiles"):
        if binning == "fisher_jenks":
            binning_method = FisherJenks(geodf[column], k=k)
        else:
            binning_method = Quantiles(geodf[column], k=k)
        bins = np.insert(binning_method.bins, 0, geodf[column].min())
        geodf = geodf.assign(__bin__=binning_method.yb)
    elif binning == "uniform":
        bins = np.linspace(
            min_value, max_value + (max_value - min_value) * 0.001, num=k + 1
        )
        geodf = geodf.assign(
            __bin__=lambda x: pd.cut(
                x[column], bins=bins, include_lowest=True, labels=False
            ).astype(np.int)
        )
    elif binning == "custom":
        if bins is None:
            raise ValueError("bins are needed for custom binning")
        bins = np.array(bins)
        geodf = geodf.assign(
            __bin__=lambda x: pd.cut(
                x[column], bins=bins, include_lowest=True, labels=False
            ).astype(np.int)
        )
    else:
        raise ValueError(
            "only fisher_jenks, quantiles and uniform binning are supported"
        )

    cmap_name = None
    norm = None
    midpoint = 0.0
    using_divergent = False
    built_palette = None

    if palette_center is not None:
        midpoint = palette_center

    if palette is None:
        if min_value < 0 and max_value > 0:
            # divergent
            cmap_name = default_divergent
            norm = MidpointNormalize(vmin=min_value, vmax=max_value, midpoint=midpoint)
            using_divergent = True
        elif min_value < 0:
            # sequential
            cmap_name = default_negative
        else:
            # sequential
            cmap_name = default_positive
    else:
        if not isinstance(palette, colors.Colormap):
            if colors.is_color_like(palette):
                if palette_type == "light":
                    built_palette = sns.light_palette(palette, n_colors=k)
                else:
                    built_palette = sns.dark_palette(palette, n_colors=k)
            else:
                built_palette = sns.color_palette(palette, n_colors=k)

    if norm is None:
        if palette_center is None:
            # norm = colors.Normalize(vmin=min_value, vmax=max_value)
            norm = colors.BoundaryNorm(bins, k)
        else:
            norm = MidpointNormalize(vmin=min_value, vmax=max_value, midpoint=midpoint)
            using_divergent = True

    if built_palette is None:
        if not using_divergent:
            built_palette = sns.color_palette(cmap_name, n_colors=k)
        else:
            middle_idx = np.where((bins[:-1] * bins[1:]) <= 0)[0][0]
            left = middle_idx
            right = k - middle_idx - 1
            if left == right:
                built_palette = sns.color_palette(cmap_name, n_colors=k)
            else:
                delta = np.abs(left - right)
                expanded_k = k + delta
                start_idx = max(middle_idx - left, right - middle_idx)
                built_palette = sns.color_palette(cmap_name, n_colors=expanded_k)[
                    start_idx : start_idx + k
                ]

    for idx, group in geodf.groupby("__bin__"):
        group.plot(
            ax=ax,
            facecolor=built_palette[idx],
            linewidth=linewidth,
            edgecolor=edgecolor,
            alpha=alpha,
            zorder=zorder,
            aspect=None,
        )

    if legend is not None:
        cbar_ax = add_ranged_color_legend(ax, bins, built_palette, **cbar_args)
    else:
        cbar_ax = None

    return ax, cbar_ax

from mapclassify import FisherJenks, Quantiles def choropleth_map( ax, geodf, column, k=10, palette=None, default_divergent="RdBu_r", default_negative="Blues_r", default_positive="Reds", palette_type="light", legend="colorbar", edgecolor="white", palette_center=None, binning="uniform", bins=None, alpha=1.0, linewidth=1, zorder=1, cbar_args={}, **kwargs, ): geodf = geodf[pd.notnull(geodf[column])].copy() min_value, max_value = geodf[column].min(), geodf[column].max() if binning in ("fisher_jenks", "quantiles"): if binning == "fisher_jenks": binning_method = FisherJenks(geodf[column], k=k) else: binning_method = Quantiles(geodf[column], k=k) bins = np.insert(binning_method.bins, 0, geodf[column].min()) geodf = geodf.assign(__bin__=binning_method.yb) elif binning == "uniform": bins = np.linspace( min_value, max_value + (max_value - min_value) * 0.001, num=k + 1 ) geodf = geodf.assign( __bin__=lambda x: pd.cut( x[column], bins=bins, include_lowest=True, labels=False ).astype(np.int) ) elif binning == "custom": if bins is None: raise ValueError("bins are needed for custom binning") bins = np.array(bins) geodf = geodf.assign( __bin__=lambda x: pd.cut( x[column], bins=bins, include_lowest=True, labels=False ).astype(np.int) ) else: raise ValueError( "only fisher_jenks, quantiles and uniform binning are supported" ) cmap_name = None norm = None midpoint = 0.0 using_divergent = False built_palette = None if palette_center is not None: midpoint = palette_center if palette is None: if min_value < 0 and max_value > 0: # divergent cmap_name = default_divergent norm = MidpointNormalize(vmin=min_value, vmax=max_value, midpoint=midpoint) using_divergent = True elif min_value < 0: # sequential cmap_name = default_negative else: # sequential cmap_name = default_positive else: if not isinstance(palette, colors.Colormap): if colors.is_color_like(palette): if palette_type == "light": built_palette = sns.light_palette(palette, n_colors=k) else: built_palette = sns.dark_palette(palette, n_colors=k) else: built_palette = sns.color_palette(palette, n_colors=k) if norm is None: if palette_center is None: # norm = colors.Normalize(vmin=min_value, vmax=max_value) norm = colors.BoundaryNorm(bins, k) else: norm = MidpointNormalize(vmin=min_value, vmax=max_value, midpoint=midpoint) using_divergent = True if built_palette is None: if not using_divergent: built_palette = sns.color_palette(cmap_name, n_colors=k) else: middle_idx = np.where((bins[:-1] * bins[1:]) <= 0)[0][0] left = middle_idx right = k - middle_idx - 1 if left == right: built_palette = sns.color_palette(cmap_name, n_colors=k) else: delta = np.abs(left - right) expanded_k = k + delta start_idx = max(middle_idx - left, right - middle_idx) built_palette = sns.color_palette(cmap_name, n_colors=expanded_k)[ start_idx : start_idx + k ] for idx, group in geodf.groupby("__bin__"): group.plot( ax=ax, facecolor=built_palette[idx], linewidth=linewidth, edgecolor=edgecolor, alpha=alpha, zorder=zorder, aspect=None, ) if legend is not None: cbar_ax = add_ranged_color_legend(ax, bins, built_palette, **cbar_args) else: cbar_ax = None return ax, cbar_ax

In [41]:

fig, ax = figure_from_geodataframe(zones, height=9, remove_axes=False)

ax, cax = choropleth_map(ax, zones.join(distancia_zonas, how='inner'), 'distancia_al_trabajo', k=5)
cax.set_title('Distancia al Trabajo [m]', loc='left')
ax.grid(True)

fig.tight_layout()

fig, ax = figure_from_geodataframe(zones, height=9, remove_axes=False) ax, cax = choropleth_map(ax, zones.join(distancia_zonas, how='inner'), 'distancia_al_trabajo', k=5) cax.set_title('Distancia al Trabajo [m]', loc='left') ax.grid(True) fig.tight_layout()

C:\Users\franc\AppData\Local\Temp\ipykernel_2608\3905406884.py:41: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  ).astype(np.int)
C:\Users\franc\AppData\Local\Temp\ipykernel_2608\1264298429.py:7: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  fig.tight_layout()

La leyenda explica que valores más claros tienen menores distancias al trabajo. Al mismo tiempo, observamos la concentración de los trabajos en la ciudad: las zonas más claras tienden a aparecer en Santiago, Providencia, Las Condes, y sus alrededores. Como consecuencia, las zonas periféricas tienen una distancia al trabajo mayor: casi toda la periferia está compuesta de zonas oscuras.

¿cómo se calculan los rangos que son representados por cada color?

Existen varias alternativas. La opción por omisión es dividir el intervalo de valores posibles en rangos uniformes (opción binning='uniform'). Otra es dividir las observaciones en cuantiles (opción binning='quantiles'). Y otra es utilizar el esquema llamado Natural Breaks, que busca la mejor manera de clasificar datos en un número específico de categorías (una especie de clustering K-Means en una dimensión). Lo podemos especificar en la función de dibujo con el parámetro binning='fisher_jenks' (se llama así porque lo inventó Jenks y lo mejoró Fisher).

Usemos la función small_multiples_from_geodataframe para comparar estas tres estrategias de codificación visual en el canal de color:

In [42]:

def small_multiples_from_geodataframe(
    geodf,
    n_variables,
    height=5,
    col_wrap=5,
    bbox=None,
    sharex=True,
    sharey=True,
    remove_axes=True,
    set_limits=True,
    flatten_axes=True,
    aspect="auto",
    basemap=None,
    basemap_interpolation="hanning",
):
    if n_variables <= 1:
        return figure_from_geodataframe(
            geodf,
            height=height,
            bbox=bbox,
            remove_axes=remove_axes,
            set_limits=set_limits,
            basemap=basemap,
            basemap_interpolation=basemap_interpolation,
        )

    if bbox is None:
        bbox = geodf.total_bounds

    # code from geopandas
    if aspect == "auto":
        if geodf.crs and geodf.crs.is_geographic:
            y_coord = np.mean([bbox[1], bbox[3]])
            aspect_ratio = 1 / np.cos(y_coord * np.pi / 180)
            # formula ported from R package sp
            # https://github.com/edzer/sp/blob/master/R/mapasp.R
        else:
            aspect_ratio = 1
    else:
        aspect_ratio = aspect

    n_columns = min(col_wrap, n_variables)
    n_rows = n_variables // n_columns
    if n_rows * n_columns < n_variables:
        n_rows += 1

    fig, axes = plt.subplots(
        n_rows,
        n_columns,
        figsize=(n_columns * height / aspect_ratio, n_rows * height),
        sharex=sharex,
        sharey=sharey,
        squeeze=False,
    )
    flattened = axes.flatten()

    if set_limits:
        for ax in flattened:
            ax.set_xlim([bbox[0], bbox[2]])
            ax.set_ylim([bbox[1], bbox[3]])

    for ax in flattened:
        ax.set_aspect(aspect_ratio)

    if remove_axes:
        for ax in flattened:
            ax.set_axis_off()
    else:
        # deactivate only unneeded axes
        for i in range(n_variables, len(axes)):
            flattened[i].set_axis_off()

    if basemap is not None:
        for ax in flattened:
            cx.add_basemap(
                ax,
                crs=geodf.crs.to_string(),
                source=basemap,
                interpolation=basemap_interpolation,
                zorder=0,
            )

    if flatten_axes:
        return fig, flattened

    return fig, axes

def small_multiples_from_geodataframe( geodf, n_variables, height=5, col_wrap=5, bbox=None, sharex=True, sharey=True, remove_axes=True, set_limits=True, flatten_axes=True, aspect="auto", basemap=None, basemap_interpolation="hanning", ): if n_variables <= 1: return figure_from_geodataframe( geodf, height=height, bbox=bbox, remove_axes=remove_axes, set_limits=set_limits, basemap=basemap, basemap_interpolation=basemap_interpolation, ) if bbox is None: bbox = geodf.total_bounds # code from geopandas if aspect == "auto": if geodf.crs and geodf.crs.is_geographic: y_coord = np.mean([bbox[1], bbox[3]]) aspect_ratio = 1 / np.cos(y_coord * np.pi / 180) # formula ported from R package sp # https://github.com/edzer/sp/blob/master/R/mapasp.R else: aspect_ratio = 1 else: aspect_ratio = aspect n_columns = min(col_wrap, n_variables) n_rows = n_variables // n_columns if n_rows * n_columns < n_variables: n_rows += 1 fig, axes = plt.subplots( n_rows, n_columns, figsize=(n_columns * height / aspect_ratio, n_rows * height), sharex=sharex, sharey=sharey, squeeze=False, ) flattened = axes.flatten() if set_limits: for ax in flattened: ax.set_xlim([bbox[0], bbox[2]]) ax.set_ylim([bbox[1], bbox[3]]) for ax in flattened: ax.set_aspect(aspect_ratio) if remove_axes: for ax in flattened: ax.set_axis_off() else: # deactivate only unneeded axes for i in range(n_variables, len(axes)): flattened[i].set_axis_off() if basemap is not None: for ax in flattened: cx.add_basemap( ax, crs=geodf.crs.to_string(), source=basemap, interpolation=basemap_interpolation, zorder=0, ) if flatten_axes: return fig, flattened return fig, axes

In [43]:

fig, axes = small_multiples_from_geodataframe(zones, 3, height=9, remove_axes=False)

for binning, ax in zip(['uniform', 'quantiles', 'fisher_jenks'], axes):
    ax.set_title(f'binning = {binning}')
    _, cax = choropleth_map(ax, zones.join(distancia_zonas, how='inner'), 'distancia_al_trabajo', k=5, binning=binning, linewidth=0.1)
    cax.set_title('Distancia al Trabajo [m]', loc='left')
    ax.grid(True)
fig.tight_layout()

fig, axes = small_multiples_from_geodataframe(zones, 3, height=9, remove_axes=False) for binning, ax in zip(['uniform', 'quantiles', 'fisher_jenks'], axes): ax.set_title(f'binning = {binning}') _, cax = choropleth_map(ax, zones.join(distancia_zonas, how='inner'), 'distancia_al_trabajo', k=5, binning=binning, linewidth=0.1) cax.set_title('Distancia al Trabajo [m]', loc='left') ax.grid(True) fig.tight_layout()

C:\Users\franc\AppData\Local\Temp\ipykernel_2608\3905406884.py:41: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  ).astype(np.int)
C:\Users\franc\AppData\Local\pypoetry\Cache\virtualenvs\db-connectors-RsEieBu8-py3.8\lib\site-packages\mapclassify\classifiers.py:1956: UserWarning: Numba not installed. Using slow pure python version.
  warnings.warn(
C:\Users\franc\AppData\Local\Temp\ipykernel_2608\3215009884.py:8: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  fig.tight_layout()

¿Cuál esquema les parece mejor (para una definición de "mejor" elegida por ustedes)? Cada uno tiene sus ventajas y desventajas.

¿Qué otro patrón observan en este mapa? ¿Qué más harían?

¿Se imaginan este mapa mezclado con los mapas anteriores? Ése es un ejercicio propuesto.

El esquema de Fisher-Jenks es utilizado con frecuencia posiblemente porque balancea los aspectos de los otros esquemas. Lo utilizaremos para la versión final del mapa, esta vez con información contextual. Además de la información territorial, agregaremos elementos típicos de los mapas cartográficos: una escala y una flecha hacia el norte.

In [45]:

grid = GeoFacetGrid(zones.join(distancia_zonas, how="inner"), height=9)

grid.add_basemap("data/scl_toner_12.tif")

grid.add_layer(
    choropleth_map,
    "distancia_al_trabajo",
    k=5,
    linewidth=0.5,
    edgecolor="black",
    binning="fisher_jenks",
    palette="RdPu",
    alpha=0.85,
    cbar_args=dict(
        label="Distancia [m]",
        height="22%",
        width="2%",
        orientation="vertical",
        location="center left",
        label_size="small",
        bbox_to_anchor=(0.0, 0.0, 0.9, 1.0),
    ),
)
grid.add_map_elements()
grid.set_title("Distancia al Trabajo")
grid.tight_layout()
plt.show()
#grid.savefig('../reports/figures/example_choropleth.png', dpi=150)

grid = GeoFacetGrid(zones.join(distancia_zonas, how="inner"), height=9) grid.add_basemap("data/scl_toner_12.tif") grid.add_layer( choropleth_map, "distancia_al_trabajo", k=5, linewidth=0.5, edgecolor="black", binning="fisher_jenks", palette="RdPu", alpha=0.85, cbar_args=dict( label="Distancia [m]", height="22%", width="2%", orientation="vertical", location="center left", label_size="small", bbox_to_anchor=(0.0, 0.0, 0.9, 1.0), ), ) grid.add_map_elements() grid.set_title("Distancia al Trabajo") grid.tight_layout() plt.show() #grid.savefig('../reports/figures/example_choropleth.png', dpi=150)

C:\Users\franc\AppData\Local\pypoetry\Cache\virtualenvs\db-connectors-RsEieBu8-py3.8\lib\site-packages\mapclassify\classifiers.py:1956: UserWarning: Numba not installed. Using slow pure python version.
  warnings.warn(
C:\Users\franc\AppData\Local\pypoetry\Cache\virtualenvs\db-connectors-RsEieBu8-py3.8\lib\site-packages\seaborn\axisgrid.py:118: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  self._figure.tight_layout(*args, **kwargs)

De esta manera el mapa nos permite caracterizar las zonas, y conjuntos de zonas, de acuerdo a la distancia que deben recorrer sus habitantes para ir al trabajo.

Queda propuesto desagregar este mapa por características de las personas, como su nivel socioeconómico.